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Abstract 

Western European coastal sand dunes are highly important for nature conservation. Communities of the creeping willow 
(Salix repens) represent one of the most characteristic and diverse vegetation types in the dunes. We report here the results 
of the first kingdom-wide fungal diversity assessment in S. repens coastal dune vegetation. We carried out massively parallel 
pyrosequencing of ITS rDNA from soil samples taken at ten sites in an extended area of joined nature reserves located along 
the North Sea coast of the Netherlands, representing habitats with varying soil pH and moisture levels. Fungal communities 
in Salix repens beds are highly diverse and we detected 1211 non-singleton fungal 97% sequence similarity OTUs after 
analyzing 688,434 ITS2 rDNA sequences. Our comparison along a north-south transect indicated strong correlation between 
soil pH and fungal community composition. The total fungal richness and the number OTUs of most fungal taxonomic 
groups negatively correlated with higher soil pH, with some exceptions. With regard to ecological groups, dark-septate 
endophytic fungi were more diverse in acidic soils, ectomycorrhizal fungi were represented by more OTUs in calcareous 
sites, while detected arbuscular mycorrhizal genera fungi showed opposing trends regarding pH. Furthermore, we detected 
numerous red listed species in our samples often from previously unknown locations, indicating that some of the fungal 
species currently considered rare may be more abundant in Dutch S. repens communities than previously thought. 



Citation: Geml J, Gravendeel B, van der Gaag KJ, Neilen M, Lammers Y, et al. (2014) The Contribution of DNA Metabarcoding to Fungal Conservation: Diversity 
Assessment, Habitat Partitioning and Mapping Red-Listed Fungi in Protected Coastal Salix repens Communities in the Netherlands. PLoS ONE 9(6): e99852. doi:10. 
1 371 /journal.pone.0099852 

Editor: Zhengguang Zhang, Nanjing Agricultural University, China 

Received March 17, 2014; Accepted May 17, 2014; Published June 17, 2014 

Copyright: © 2014 Geml et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits 
unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 

Data Availability: The authors confirm that all data underlying the findings are fully available without restriction. Data are available from the European 
Nucleotide Archive under the accession number: ERP001713. 

Funding: This research was sponsored by a Next-Generation Sequencing grant from the Naturalis Biodiversity Center to JG and MEN and by the Kits van Waveren 
Postdoctoral Fellowship to JG. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. 

Competing Interests: The authors have declared that no competing interests exist. 

* Email: Jozsef.Geml@naturalis.nl 



Introduction 

Fungi represent one of the largest groups of living organisms. 
However, fungi are still poorly understood and appreciated 
compared to plants and animals and our knowledge of fungal 
diversity lags far behind [1]. Approximately 100,000 species of 
Eumycota (the true Fungi) have been described. Already well 
before the routine use of DNA sequencing in fungal biodiversity 
assessments, their true diversity was estimated to be around 1.5 
million species [2] , while more recent estimates suggests there may 
be 0.7 to 5 million species [3], [4]. Despite the differences among 
various estimates, it is clear that we currently know only a fraction 
of the total fungal biodiversity. 

The task of discovering a significant portion of yet unknown 
species has only recently become possible with the advent of high- 
throughput DNA sequencing of environmental samples that has 
enormous potential to further boost data acquisition in biodiversity 
research [5], [6], [7], [8], [9], [10], [1 1], [12], [13], [14], [15], [16]. This 



methodology involves identification of multiple species from 
environmental samples using targeted loci specifically selected 
with the purpose of identification. These loci are typically the same 
ones that are used in specimen-based large-scale DNA barcoding 
efforts (e.g., [17], [18]), hence providing both appropriate power of 
resolution for the group of interest and maximum amount of 
reference sequences available from vouchered specimens. 

The fungal diversity assessment presented here focuses on Salix 
repens L. sand dune communities along the North Sea coast, 
because these areas are highly important for nature conservation, 
water resource management, and recreational purposes. Coastal 
sand dunes include an outstanding ecological diversity in terms of 
environmental heterogeneity and species composition [19]. In the 
Netherlands, they are home to more than half of all higher plant 
species, in spite of covering only 1% of the land area [20], [21]. In 
Europe, coastal dune areas have been reduced by ca. 70% in the 
last century [22], and biodiversity in the remaining dune 
ecosystems is highly endangered [23]. In the Netherlands, 
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however, coastal dunes have been relatively well preserved, with 
less than 10% area loss, largely due to their roles in drinking water 
management and as barriers against the sea [21]. These dunes 
represent one of the last nutrient-poor habitats in the Netherlands 
and most plant species featured on the Dutch Red List now have 
their main distribution in these coastal areas. Moreover, the 
reserves featured in this manuscript are part of Natura 2000, an 
ecological network of protected areas in the territory of the 
European Union (http://www.natura2000.nl/pages/kaartpagina. 
aspx). 

The creeping willow (Salix repens) is a widespread shrub in 
western Europe and inhabits a wide range of vegetation types. It 
forms both ectomycorrhiza and arbuscular mycorrhiza [2 4], [25], 
often supporting a diverse community of symbiotic fungi in its 
extensive root system [26]. It is a particularly characteristic 
component of coastal sand dune ecosystems where it is often the 
only ectomycorrhizal (ECM) host plant [27] , [26] . Therefore, it has 
been repeatedly noted that the presence of S. repens adds wealth to 
the overall fungal diversity in sand dune ecosystems [25], [26]. 
Somewhat surprisingly, species composition of fungi in S. repens 
communities have not been previously investigated using molec- 
ular methods to our knowledge. To date, characterizations of fungi 
in S. repens assemblages have been based on morphology-based 
identification of fruitbodies and ECM root structures. For 
example, in one of the early studies comparing above- and 
below-ground diversity of ECM fungi, Van der Heijden et al. [25] 
found fruitbodies of 78 ECM taxa belonging to 12 genera and 
identified fifteen ECM root morphotypes in S. repens communities 
in the Netherlands. Watling [26] reported almost 80 species of 
macrofungi from S. repens beds of the northern Scottish islands 
based on fruiting records. Despite these outstanding pioneer 
works, important gaps remain in our knowledge about total fungal 
diversity in this important ecosystem. 

In this study, we collected soil samples at multiple locations in 
strictiy protected coastal areas and conducted a large-scale 
biodiversity assessment of soil fungi inhabiting the rhizosphere of 
S. repens using pyrosequencing techniques. The aims of this work 
were 1) to characterize fungal communities found in S. repens 
vegetation; 2) compare the effects of environmental factors (e.g., 
soil pH, geological district) on communities through multiple sites; 
and 3) to supplement sporocarp-based mapping data points of rare 
fungi with soil-based data points. Our results provide the first 
kingdom-wide fungal diversity assessment in coastal sand dune 
communities and demonstrate the potential use of high-through- 
put DNA sequencing of environmental samples to fungal 
conservation. 

Materials and Methods 

The Study Region 

The following environmental agencies granted permission for 
the fieldwork (with sampled field sites): Amsterdamse Water- 
leidingduinen (Haarlem-Bentveld), Dunea (Wassenaar-Meijendel), 
Natuurmonumenten (Goeree, Callantsoog-Zwanenwater), and 
Staatsbosbeheer (Terschelling). GPS coordinates of the sampling 
sites are listed in Table 1. The field sampling only included soil 
and did not involve endangered or protected species. The selected 
sites represented sand dune S. repens beds (European Environment 
Agency habitat code: 2170) in an extended area of joined nature 
reserves located along the North Sea coast of the Netherlands, 
representing a south-north transect (Figure 1, Table 1). This 
coastal dune system is characterized by an oceanic climate with 
mean annual temperature of 9-10°C (July, 1 6-1 7 °C; January, 2- 
3°C) and precipitation of 675-800 mm [21]. Sand-lime content 



varies between 0.1 and 10% [21], with generally calcareous soils in 
the southern Renodunal district (represented by the sampling sites 
Goeree, Meijendel, and Haarlem-Bentveld) and more acidic soils 
in the northern Wadden district (Zwanenwater, Terschelling). In 
most areas, the height of the dunes reaches 15-25 m [28]. Current 
management policies are directed toward multiple purposes, such 
as drinking water supply, defense against storm floods, recreation, 
and nature conservation [29]. Soil samples were taken at sites 
monitored by Natuurmonumenten (Goeree and Zwanenwater), 
Dunea (Meijendel), Staatsbosbeheer (Terschelling), and the 
Amsterdamse Waterleidingduinen (Haarlem-Bentveld). In March 
of 2010, 20 soil cores, 2 cm in diameter and ca. 20 cm deep, were 
taken at each site in a way that cores were at least 2 m from each 
other to minimize the probability of sampling the same genet 
repeatedly. For every plot, the 20 cores were pooled, resulting in a 
composite soil sample for each plot. Soils were deposited in 50 mL 
Falcon tubes and stored at -80°C until lyophilization. 

DNA Extraction, PCR, and Pyrosequencing 

We follow the recently proposed guidelines for standardizing the 
description of next-generation sequencing datasets in publications 
as proposed by Nilsson et al. [30]. Genomic DNA was extracted 
from 8 g of lyophilized soil from each composite sample using the 
Mo Bio Powersoil kit following the manufacturer's instructions. 
We used high-throughput tag-encoded 454 GS-FLX amplicon 
sequencing to assess fungal diversity in the selected plant 
communities. Amplicon libraries were performed using a combi- 
nation of tagged primers designed for the variable ITS region, as 
recommended for the tag-encoded pyrosequencing method [31]. 
Genomic DNA samples were amplified using the fungal-specific 
primer pair ITS1F (5'-ACTTGGTCATTTAGAGGAAGTAA) 
and ITS4 5'-) \uuT( :( !T( ( GCTTATTGATATGC . where X 
and Y represent the two pyrosequencing primers (CCTATCCC- 
CTGTGTGCCTTGGCAGTCTCAGT and CCATCTCATC- 
CCTGCGTGTCTCCGACTCAGA) and xxxxx represents the 
barcodes designed for sample identification. Barcodes differed 
from all other barcodes by at least two nucleotides. Four PCR 
reactions were performed per sample using the following 
temperature program: 95°C/2 min, 25 cycles of 95°C/30 s, 
54°C/1 min, 72°C/2 min, and 72°C/10 min. One ul of DNA 
template was used for the 40 JLXl PCR reaction containing 25.6 ul 
of MQ, water, 4 ul of lOx buffer, 1.5 ul dNTP's (2.5 mM), 1.5 ul 
of reverse and forward primers (10 mM), 4 ul MgCl 2 (50 mM), 
0.5 ul BSA (10 mg/ml) and 0.4 ul BIOTAQ, polymerase (5 U/ul). 
PCR products were pooled for each sample, purified and 
sequencing-adapters/tags were integrated via MegaPrimers 
(PCR). The amplicon length and concentration were estimated, 
and an equimolar mix of all amplicon libraries was used for 
pyrosequencing (with the ITS4 primer). Tagged samples were 
pooled and sequenced using the Genome Sequencer GS FLX 454 
System (454 Life Sciences/Roche Applied Biosystems, Nutley, NJ, 
USA) at GATC Biotech AG (Konstanz, Germany) and the Leiden 
Genome Technology Center (Leiden, the Netherlands). 

Bioinformatic Work 

Pyrosequencing resulted in 688 434 reads following quality 
checks and trimming, with median read length of 284 bp. The raw 
data have been submitted to the European Nucleotide Archive 
(ERP001713). We used a parallel version of MOTHUR v. 1.32.1 
[32] installed at the University of Alaska Life Sciences Informatics 
Portal for sequence analyses. We filtered sequences by quality 
score with a sliding window approach as in Brown et al. [33]. The 
FASTQ, files were converted to FASTA and QUAL files, and the 
sequences were subjected to quality filtering whereby each 
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Figure 1. Sampling sites along the coast of the Netherlands. Sites corresponding to the Renodunal and Wadden geological districts are 
marked with black and white squares, respectively. 
doi:1 0.1 371 /journal.pone.0099852.g001 



sequence was screened for thresholds for average Phred score of 
Q&25 in a sliding window of 50 bp (qwindowaverage = 25, 
qwindowsize = 50), no ambiguous bases (maxambig = 0) and 
homopolymers no longer than 8 bp (maxhomop = 8). Sequences 
shorter than 150 bp or longer than 400 bp were omitted from 
further analyses (mintlength= 150, maxlength = 400). A total of 
369 181 sequences remained after quality filtering and trimming. 
Because next-generation sequencing libraries generally vary in 
size, we normalized the number of sequences for all samples prior 
to OTU clustering, as recommended by Gihring et al. [34], to 
ensure that estimators across all samples are comparable. For this 
purpose, we randomly subsampled the number of trimmed and 
quality-filtered reads to the size of the smallest library (30 497). 



The resulting 304 970 sequences served as input for operational 
taxonomic unit (OTU) clustering using OTUPIPE [35] based on 
97% sequence similarity. Such 97% sequence similarity cut-off 
value is routinely used in fungal community studies as a proxy for 
species, accounting for intraspecific variation and errors generated 
by PCR (e.g., [36], [37], [38], [7], [39], [40], [9], [12], [15], [41]. We 
simultaneously removed 13 279 putatively chimeric sequences, 
using a curated dataset of fungal ITS sequences as reference 
dataset [42]. Pyrosequencing datasets comprising orders of 
magnitude more sequences per taxa than traditional methods 
tend to contain a higher proportion of artifactual OTUs, because 
sequencing errors tend to be unique, while the accumulation of 
real taxa starts to level off [43], [12]. Because most pyrosequencing 
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Table 1. Sampling sites of coastal Salix repens communities used in this study. 



Site 


Code 


Geological district, habitat 


Latitude, longitude 


PH 


OTUs 


Goeree: Middelduinen 


R-G1 


Renodunal, wet 


51.827238; 3.943405 


6.01 


971 


Goeree: Westduinen 


R-G2 


Renodunal, wet 


51.810850; 3.861609 


6.73 


934 


Wassenaar: Meijendel 1 


R-M1 


Renodunal, dry 


52.137439; 4.341960 


7.3 


691 


Wassenaar: Meijendel 2 


R-M2 


Renodunal, dry 


52.142707; 4.333034 


7.42 


790 


Haarlem-Bentveld: Amsterdamse Waterleidingduinen 1 


R-H1 


Renodunal, wet 


52.338001; 4.531978 


7.11 


727 


Haarlem-Bentveld: Amsterdamse Waterleidingduinen 2 


R-H2 


Renodunal, dry 


52.337312; 4.533910 


6.56 


660 


Callantsoog: Zwanenwater 1 


W-Z1 


Wadden, dry 


52.822073; 4.709229 


5.71 


867 


Callantsoog: Zwanenwater 2 


W-Z2 


Wadden, wet 


52.821516; 4.708264 


5.32 


762 


Terschelling: Groene Strand 


W-T1 


Wadden, wet 


53.364741; 5.201769 


6.39 


718 


Terschelling: West aan Zee 


W-T2 


Wadden, dry 


53.398326; 5.269489 


5.78 


804 



Code, geological district, habitat type (wet dune slack vs. dry dune side and top), geographic coordinates, soil pH, number of quality-filtered sequences and the number 
of 97% ITS2 sequence similarity OTUs are displayed for each site. 
doi:1 0.1 371 /journal.pone.0099852.t001 



singletons tend to be artifactual and such sequencing errors can 
overestimate the diversity of 'rare taxa' (e.g., [44], [12]), we opted 
to be conservative and excluded all singletons from further 
analyses. For identification, sequences were subjected to 
USEARCH [45] against the latest release of quality-checked 
UNITE+INSD fungal ITS sequence database containing both 
identified and unidentified sequences, many of which are assigned 
to Species Hypothesis groups as defined by experts using 
phylogenetic evaluations [46]. OTUs that did not have at least 
80% similarity over at least 100 bp to any fungal sequence were 
excluded from further analyses. Taxonomic affiliations were based 
on the current Index Fungorum classification as implemented in 
UNITE. 

Comparing Fungal Communities among Sampling Sites 

We used the 97% ITS sequence similarity OTUs inferred 
earlier as input data for the various ordination methods. To 
visualize the variation in the occurrence of OTUs across sites, we 
carried out ordination using non-metric multidimensional scaling 
(NMDS) with the quantitative version of Sorenson similarity (Bray- 
Curtis index) with 500 iterations in PC-Ord [47] . Read abundance 
in 454 sequencing has been shown to be approximately 
quantitative only within species, while between-species compari- 
sons can be biased by innate sequence structure [39]. Because of 
such uncertainties regarding the reliability of read abundance, we 
carried out the ordination analyses with square-root transformed 
abundance to moderate the influence of OTUs with high sequence 
counts, while maintaining some approximation of abundance that 
often reflects ecological significance. The solution with the lowest 
stress was derived from 250 runs using real data and was then 
subjected to 250 randomized runs using a Monte Carlo test to 
evaluate the probability of final the NMDS patterns being greater 
than chance occurrences. Orthogonal rotation of the resulting 
NMDS solution was used to maximize correlation between the 
strongest environmental variables (as indicated by Pearson r 
values) and major axes. The solution with the lowest number of 
dimensions was selected when the decrease in the final stress was 
<5 by adding another dimension [47]. The Pearson correlation 
coefficient (r) values between environmental and fungal commu- 
nity variables and axes 1 and 2 were calculated. Soil pH and the 
number of OTUs representing various taxonomic groups per site 
were included in a site environmental factor matrix as quantitative 
variables. For this latter, we used only taxonomic orders and 



genera represented by at least 3 OTUs were included for a site as 
quantitative variables (see Table 2 for full list). Geological district 
(Renodunal = 1, Wadden = 0) was included as categorical variable. 
Following NMDS ordination of sites, we examined the Pearson 
correlation values between the community ordination axes and 
these variables. We also tested whether fungal communities were 
statistically different across the sites using a multiresponse 
permutation procedure (MRPP), also in PC-Ord. For the MRPP 
analyses, three categorical variables were tested: geological district, 
as described above, habitat type (wet dune slack vs. dry dune side 
and top), and soil pH (acidic = 0, non-acidic = 1 , based on a 
pH = 6.5 threshold following the USDA soil classification for 
acidic and non-acidic soils). Finally, we determined any prefer- 
ences of individual OTUs for pH categories using indicator species 
analyses. 

Results 

OTU Diversity and Taxonomic Affiliations 

Out of the 688 434 original reads, 402 500 sequences passed the 
various filtering steps. After normalizing the library size across all 
samples, 304 970 sequences were assembled into 7156 OTUs of 
which 1788 were singletons. After excluding all singletons and 
OTUs with <80% similarity or < 1 00 bp pairwise alignment 
length to a fungal sequence as performed by USEARCH, 3611 
OTUs were retained for subsequent analyses. The number of 
OTUs at a given site ranged from 660 (Amsterdamse Water- 
leidingduinen, R-H2) to 971 (Goeree: Middelduinen, R-Gl) 
(Table 1). 

Ascomycota was the dominant phylum and accounted for 
33.18% of the OTUs, followed by Basidiomycota (22.73%) 
Glomeromycota (5.29%), Mucoromycotina, incertae sedis (1.94%), 
and Chytridiomycota (0.53%), while unidentified fungal OTUs 
with most similar sequences to other environmental sequences 
without assignment to phylum accounted for 36.33% (Figure 2a). 
In Ascomycota, the ranking of taxonomic orders based on the 
number of representative OTUs was as follows: Helotiales (171), 
Pleosporales (136), Hypocreales (114), Pezizales (92), Archaeorhi- 
zomycetales (71), Chaetothyriales (64), Capnodiales (51), Sordar- 
iales (45), Eurotiales (34), Xylariales (19), Diaporthales (17), 
Verrucariales (13), Lecanorales (13), Thelebolales (12), Botryopsh- 
paeriales (10), followed by numerous other orders with <10 
OTUs. In addition, there were 49 ascomycete OTUs with incertae 
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(b) 
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Figure 2. Taxonomic groupings of fungi retrieved from the samples. Proportional distribution of the 1022 unique BLASTN hits among 
different phyla and taxonomic orders in Basidiomycota and Ascomycota across all sampling sites based on sequences most closely related to 97% ITS 
OTUs in our sample. 
doi:1 0.1 371 /journal.pone.0099852.g002 



sedis classification, representing mostly mitosporic, asexual fungi 
(Figure 2b). In Basidiomycota, Agaricales was by far the most 
species-rich order in the sample with 335 OTUs, followed by 
Thelephorales (116), Cantharellales (69), Sebacinales (67), Tre- 
mellales (30), Sporidiobolales (26), Polyporales (19), Filobasidiales 
(16), Russulales (15), Trechisporales (11) and numerous other 
orders with <10 OTUs each, while there were only 2 
basidiomycete OTUs with uncertain taxonomic placement (i.e, 
incertae sedis) (Figure 2c). 

There were 354 OTUs with >97% sequence similarity to 
publicly available sequences with species identification (Table SI). 
Although most of these represented microfungi, we also detected 
many species that form macroscopic fruitbodies (e.g., mushrooms, 
bracket fungi, cup fungi, earth tongues, truffles etc.), among which 
we found numerous red listed species in our samples, often from 
previously unknown locations, based on mapping data and red list 
categories obtained from the Fungal Mapping website of the 
Dutch Mycological Society (http://www.verspreidingsatlas.nl/ 
paddenstoelen). 

Comparing Fungal Communities among Sampling Sites 

NMDS analyses resulted in a 2-dimensional solution with a final 
stress of 0.03062 and a final instability <0.00001. The Monte 
Carlo test results indicated that this 2-dimensional solution using 
the real data was significantly better than chance occurrences 
[p = 0.008). The two axes explained the majority of variability in 
the sampled fungal communities (axis 1: r 2 = 0.855; axis 2: 
r 2 = 0.049; total f = 0.904; orthogonality = 79.5%). The NMDS 
ordination plot was orthogonally rotated by pH to visualize 
correlations between pH and fungal community composition in 
general and the taxonomic groups in particular. The Pearson 
correlation coefficient (r) values between all fungal community 
variables and axes 1 and 2 are displayed in Table 2. As in the 
fungal community ecology paper of Rogers et al. [48], we 
considered variables with \r\ SO. 5 values for axis 1 important 
for characterizing changes in fungal community structure along a 
pH gradient and these variables were superimposed on the 
ordination plot as direction and strength vectors (Figure 3). 

There was a strong structuring of the occurrence of fungal 
OTUs according to soil pH, that co-correlated with the geological 
district. Sites belonging to the Wadden district were clustered on 
the left side of the ordination plot with strongly negative axis 1 
values. Even though the Renodunal district sites were more 
scattered, the Wadden and Renodunal sites formed non-overlap- 
ping clusters (Figure 3). However, one Renodunal site (Goeree: 
Middelduinen, R-Gl), that had an acidic soil pH value 
uncharacteristic of the Renodunal sites, was appeared to be more 
similar to the Wadden communities than to other Renodunal sites. 
MRPP tests confirmed the importance of soil pH in shaping fungal 
community composition (effect size A = 0.04259, probability 
p = 0.00463), while the effect of geological district was somewhat 
weaker (A = 0.02786, p = 0.03 1 57). Habitat type (wet dune slack vs. 
dry dune side and top) did not seem to have a significant influence 
on community composition (A = 0.01 1 13, p — 0.17961). 

The joint plots illustrate the strength and direction of correlation 
of variables to each ordination axis (Figure 3). The correlation of 
soil pH with axis 1 was strong (r= 0.7 1 8) . At the phylum level, 
only Chytridiomycota showed correlation with axisl (r = —0.585) 
and was more diverse in the acidic sites, while diversity in other 



phyla did not show significant trends. At the taxonomic order 
level, the ascomycete Archaeorhizomycetales, Helotiales, Euroti- 
ales, Saccharomycetales, and Verrucariales and the basidiomycete 
Trichosporonales, Auriculariales, Trechisporales were found to be 
more diverse in the acidic sites (Figure 3, Table 2). Orders showing 
greater richness in the alkaline soil samples included the 
ascomycete Pezizales, Sordariales, Lecanorales, Pleosporales, 
Onygenales and the basidiomycete Thelephorales. At the genus 
level, the diversity in ascomycete genera Archaeorhizomyces, Lachnum, 
Penkillium, and Verticillium and basidiomycete genera Trichosporon, 
Thekphora, Clavulinopsis, Cortinarius, Clavulina, and Mycena were 
negatively correlated with pH, while the ascomycete Phoma, 
Exophiala, Podospora, Cladonia, Teratospora, Geopora, Cadophora, Alter- 
naria, Scolecobasidium and the basidiomycete Tomentella, Inocybe, 
Panaeolus, and Tanatephorus were more diverse in the alkaline 
samples (Figure 3, Table 2). 

There were 96 OTUs identified as significant (p<0.05) 
indicators for a certain pH category. Of these, 82 were indicators 
for acidic and 14 for alkaline sites. Several indicators belonged to 
root-associated fungi. For example, an unidentified root-endo- 
phytic Phialocephala and the ECM Russula pascua (SHI 17152. 05FU), 
Thekphora albomarginata, Thekphora sp. (SHI 12629. 05FU), and 
Tulasnella sp. (SHI 17022. 05FU) were associated with acidic sites, 
while the ECM Inocybe dunensis (SHI 107 18.05FU) and an 
unidentified Tomenklla were characteristic of the alkaline sites. In 
addition, there were several saprobic that were found almost 
exclusively in the acidic soil samples, such as Cryptococcus sp. 
(SH117277.05FU), Fusarium sp. (SH106152.05FU), various Mor- 
tkrella species, Mycena simian (SHI 02569. 05FU), Penkillium adamet- 
zii, Rkkenella fibula (SHI 09302. 05FU), Schizoblastosporion starkeyihen- 
ricii, Umbelopsis ramanniana (SH100505.05FU), and a Verticillium sp. 
(SHI 1 1596. 05FU), some of which can be potentially pathogenic. 
An unidentified species of the nematode pathogen Pochonia genus 
(SH106625.05FU) also showed strong preference for acidic sites. 
There were only two saprobic species identified as indicators for 
the alkaline soils: Fusarium torulosum and a Podospora species. The full 
list of indicator OTUs are shown in Table 3. 

Discussion 

Acidic vs. Non-acidic Sand Dune Fungal Communities 

While environmental factors influencing the distribution of 
plants and animals have been thoroughly studied, environmental 
factors controlling the spatial patterns and community composi- 
tion of soil microorganisms are still poorly understood [49]. 
Recent molecular studies on various microbial groups (mosdy 
bacteria) have started to explore the distributional patterns of soil 
communities, however, many of them applied DNA-fingerprinting 
techniques or phospholipid analyses that do not permit the 
taxonomic surveys and comparisons of different communities 
[50], [51]. Soil pH is generally implicitiy considered important for 
shaping fungal communities with some supporting evidence from 
experimental studies mosdy with AM fungi [52], [53] and, in 
recent years, from joint pyrosequencing analyses of bacterial and 
fungal communities in agricultural soils and land use types 
[54], [49]. Unfortunately, both of these latter sequencing studies 
relied on nuclear small-subunit (SSU) ribosomal genes, meaning 
that their power of taxonomic resolution was at family-level groups 
at best. Therefore, information on the effect of soil pH on the 
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Figure 3. Non-metric multidimensional scaling (NMDS) ordination plot for fungal communities from Salix repens sites. Ordinations 
were based on the square-root transformed abundance of 97% ITS sequence similarity OTUs. Labels, localities and descriptions of the sampling sites 
are given in Table 1. Vectors are shown for soil pH, total OTU richness, and for the most diverse taxonomic orders (A) and genera (B) that correlated 
with the ordination axes at /^>0.200. 
doi:1 0.1 371 /journal.pone.0099852.g003 



distribution of specific fungal taxa and on the composition of their 
communities in natural habitats remains surprisingly scarce. 

Our large-scale comparison of various acidic and non-acidic 
coastal dune Salix repens beds clearly indicates that soil pH strongly 
correlates with the composition of fungal communities. In NMDS 
ordinations, soil pH had high Pearson correlation values with the 
dominant ordination axis and overrode the effect of the geological 
district (Wadden vs. Renodunal) on community composition in 
both the NMDS and the MRPP analyses (Figures 3 and 4). This is 
in agreement with the findings of Van der Heijden et al. [25] that 
soil pH significandy contributed to the variation explained in 
ECM morphotype composition of communities associated with S. 
repens. In our NMDS ordinations and MRPP analyses, Inocybe and 
Tomentella showed strong preference for alkaline soils, while 
Clavulina, Cortinarius, Thelephora, and to some extent Laccaria, were 
more diverse in acidic sites (Figure 3, Table 2). It is worth noting, 
that because our study focuses on natural systems, it is impossible 
to determine based on our current data whether the communities 
are structured directly by pH or indirectly via the interaction of 
pH with soil communities and various environmental factors (e.g., 
nutrient availability, organic C characteristics, vegetation differ- 
ences etc.). Because results from previous studies using a pH 
gradient with controlled variables suggest no or only weak effect of 
pH on community structure [54], [49] and because many fungal 
species have a relatively wide pH optimum (e.g. [55], [56]), it is 
likely that the observed correlation of pH with community 
composition is mainly indirect, e.g., via nutrient availability and 
altered competitive interactions between soil fungi and bacteria 
[57], and other soil biota. While our data suggest that the effect of 
soil pH on fungal community composition may be stronger in 
relatively undisturbed sites than in agricultural settings, disentan- 
gling causal relationships is beyond the scope of our work. 

Root-associated Fungi 

Although taxa detected by our sequencing efforts featured 
groups from a wide variety of ecological guilds, including many 
saprobic and/or plant pathogenic fungi (e.g., Fusarium cereale, 
Gibberella zeae, Phoma multirostrata) as well as animal or fungal 
pathogens (e.g., Hypomyces chrysospermus, Metarhizium anisopliae, M. 
jlavoviride and other Hypocreales), we focus our discussion on 
groups that are most relevant to current conservation priorities, 
such as fungi with macroscopic fruitbodies and/ or with symbiotic 
associations with S. repens or with red listed plants (e.g. orchids) in 
these habitats. 

Previous studies based on sporocarps and ECM root structures 
have shown that Salix repens communities generally harbour 
numerous ectomycorrhizal basidiomycetes [25], [26]. This was 
confirmed by our large-scale soil sequence data, as we found 418 
ECM OTUs. The majority of these showed similarity to 
basidiomycete genera known to form associations with willows 
(e.g. [58], [59], [15]), such as Inocybe (95 OTUs), Thelephora/ 
Tomentella (83), Cortinarius (36), Hebeloma (10), Russula (7), Laccaria 
(7), Clavulina (6), Lactarius (5), Paxillus (2). Also, there were 4 OTUs 
[Amanita pantherina, Hydnum rufescens, Suillus luteus and Tulasnella sp.) 
that represent species that are known to form ECM associations 
with other host trees found near the sampling sites, such as oak, 
pine or poplar. In addition, we found 52 OTUs belonging to the 



Sebacinaceae, a taxonomically difficult group with diverse 
ecological roles, including various forms of mycorrhizae (ECM, 
orchid, ericoid, arbutoid) and other associations with plant roots 
[60]. We also identified several ectomycorrhizal ascomycetes: 
Cenococcum geophilum, one of the most widespread and abundant 
ECM fungus with a wide host range [59], Geopora arenicola and G. 
cervina, a genus known to occur in sand dunes [61], Helvetia 
compressa and H. maculata, and an unidentified Trichophaea species. 
In addition, we found four truffles that are likely symbionts of 
Quercus that is widespread in the dunes: Tuber anniae, and three 
unidentified Tuber species. 

Salix repens is one of the very few plants known to form both 
ECM and AM [2 4], [62]. In dry coastal foredunes, where 
phosphorus is the main limiting nutrient, AM fungi and their 
hosts generally predominate, while nitrogen is the main limiting 
element in dune slacks, resulting in the dominance of ECM taxa 
[62]. With the ability to form both ECM and AM associations, 
Salix repens is particularly well-suited to survive in the highly 
dynamic dune ecosystems and is able to colonize habitats that 
range from dry to wet and from acidic to calcareous [62]. We 
detected 191 glomeromycete OTUs in our samples, representing 
Glomus (53), Rhizophagus (10), Archaeospora (6), Acaulospora (4), 
Claroideoglomus (4), Scutellospora (4), Cetraspora (3), Paraglomus (3), 
Entrophospora (2), Diversispora (1), Racocetra (1), and 100 unidentified 
glomeromycete OTUs. Because the majority of land plants are 
known to form AM, it is likely that the majority of these taxa are 
also associated with a wide range of plants beside Salix repens. 

We also found numerous OTUs with very high sequence 
similarity to known dark septate root endophytes (DSE) (e.g., 
Cadophora, Cladophialophora, Exophiala, Leptodontidium, Meliniomyces, 
Phialocephala, Phialophora spp.). Although DSE colonization and 
diversity have previously been shown to be higher in acidic 
conditions in general [63], [64], our data suggest that different 
genera may prefer different pH conditions. For example, Cadophora 
and Exophiala were represented by more OTUs in the alkaline 
sites, while Meliniomyces was somewhat more diverse in samples 
with lower pH (Table 2). 

The sampling areas included in our study harbor some of the 
richest orchid populations in the Netherlands. The diversity of 
orchid mycorrhizae may be particularly relevant for nature 
conservation as they provide symbioses to this unique group of 
plants with high public appeal and conservation value. The 
majority of orchid mycorrhizal symbionts belong to the Cerato- 
basidiaceae, Sebacinaceae and Tulasnellaceae [65], [66], although 
ECM fungi have also been reported to be associated with orchids, 
albeit mainly with non-photosynthetic, mycoheterotrophic species 
that often use the ECM fungus as a channel to obtain 
carbohydrates from the ECM host tree [65], [6 7]. Among the first 
group of orchid symbionts that are typically associated with 
photosynthetic orchids, we found 21 Ceratobasidiaceae and 12 
Tulasnellaceae OTUs in addition to the 52 Sebacinaceae OTUs 
mentioned above. Among other ECM lineages, the basidiomycete 
Inocybe, Tomentella, and the ascomycete Tuber species have been 
found to form mycorrhizal associations with the mycohetero- 
trophic orchid Epipactis helleborine [6 7], [68]. This orchid is 
represented by two subspecies in the Netherlands, among these, 
E. h. ssp. helleborine is widespread in a variety of habitats, while E. h. 



PLOS ONE | www.plosone.org 



10 



June 2014 | Volume 9 | Issue 6 | e99852 



Diversity and Ecology of North Sea Dune Fungi 



■d 

6 



i " « 



— <* £ f £ 



y 4-< *j *-< +j *j +j 



E E 



5EEEEEEEE5 



E E 



fU (D TO '+j 

-f -f o o 



TOTO ^ fo fo fo fo fo 



E E E 



E £ E E E 



£ £ £ E 



q V u u u 



£ E 



^EEEEEEE 



ri. 












Q. 


Gi- 


ti- 


ti. 


to 


rc 


re 


ro 


re 


re 


o 


O 


o 


O 


o 


o 


u 
>N 


u 
>. 


u 
>. 


U 

>, 


u 
>. 


u 
>, 


E 


E 


E 


E 


E 


E 


o 


o 


o 


o 


o 


o 


u 


u 


u 


u 


u 


u 


< 


< 


< 


< 


< 


< 



P o o o o o 



-1= ^ 



I I ° 



E E 



§ § § § 



m m u-i lo lo lo 



<— co co yo \D 



rM i— N 



CO CO ^ 



X 

in 



.-00«-«3«-«-i-i- 



i/ivii/ii/it/ii/it/iui 



.Q 



O i- i- CO t- 



CO CO t— o 



<N m <N ro 



m m rsi m 



(N (N cn m 



<- <- 



C0 vD 0> 



Cft (Ji ^ d 



m m co o\ <j\ 



CO CT\ a\ CTi 



Oi i- Ol i- 



CO CO CO 



a 

X 



n o> ,- 



O CO CO 



C7\ CJ\ Q\ 



ro m J ^ 



0> 0> lo 



"1 r- £ 



r>i m t— t— 



m m in io io cc 

!> is <> *- ^ 



a 
x 



o a s 



o 



VO CO CO 



s s 



CO VO CO CO CO 

o «- o o 



.E u 



X 



D 

o 



<— i- ro m 



vo m m \0 



LTI v£> \£> \D 
Ol CO t 



' — < — vD 
CO IN T fN t 
r*. «— co co i— 



PLOS ONE | www.plosone.org 



11 



June 2014 | Volume 9 | Issue 6 | e99852 



Diversity and Ecology of North Sea Dune Fungi 



c 
o 
u 

m 

ii 
.q 
<c 
\- 



E E E 



EEEEEEEEE 



E E E E 



E E E E 



CO CO CO CO 



E E 



■5- -S *■> +* *■> 4-* 4-* +- 1 



•O SI Si 



^ u u i: ai 
o o o o 
to in m in 



u u u u 



rM «- 

LT) r~~- v£) ^ 

>— rN <— <— 



t- o 



\o en \o 



cs 



LnorMrNOi-^or-j 
r\ir\ir\irNrNmrNr\| 



\D CO 0^ 



rM ro cn co 



Oi ifl CO Oi CO 
CT\ 0*i 0*t 0*i 0*i 



0*i CO vO 



Lovo^i^yDcoo^r-^ 
o^o^o^o^o^o^o^o^ 



C( iO o 



CTi CTi CO <J\ 



ro fN r- 



10 \o ^ iO to 
t- r-^ t- m ^ 



(N <— ,_ 
~ £ £ 



J, co oo in rs 



a 



u_ u_ u_ O 



CO vO vD 0*1 O"* 



a 

X 



o> m \o 



LD u-i co m 



CO Oi CO 



u .E u 



>— m 



a 9. Q 



O r- r- 



i— m m 



i— CO i— 



v£> CO Ci Ci 



ro ro ro ro 



PLOS ONE I www.plosone.org 



12 



June 2014 | Volume 9 | Issue 6 | e99852 



Diversity and 



Ecology of North Sea Dune Fungi 



c 
o 
u 



m m <n m 



r-i r\| r\| eg 



VD LT) i — 00 



•- CT^ 



CTi IT) lo CFi 



en 



rsi m o> ro 






=> 
13 



CO CO m oo rM 



a 

Q 



oo ^ 



=^ =; < ^ ^ 



a 



a 

D 



\o oo o 



CO 00 00 CO 



I 3 

a. ™ 

T 00 
QJ O 

u q 



on o 
aj Q. 



PLOS ONE | www.plosone.org 



13 



June 2014 | Volume 9 | Issue 6 | e99852 



Diversity and Ecology of North Sea Dune Fungi 




Figure 4. Examples for OTUs found in our soil samples showing high sequence similarity with rare, red-listed fungi. Sporocarp records 
(before and after 1990) are based on the mapping database of the Dutch Mycological Society (http://www.verspreidingsatlas.nl/paddenstoelen). 
doi:1 0.1 371 /journal.pone.0099852.g004 
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ssp. neerlandica is restricted to the dunes and is endangered [68]. 
This latter rare subspecies is present in most of our sampling sites 
and we found multiple species of all three ECM genera reportedly 
associated with E. helkborine. Tomentella and Inocybe were the most 
species-rich genera in our samples and are important symbionts 
for Salix repens as well, as specified earlier. Because E. h. ssp. 
neerlandica is mostly found near Salix repens [68], [69], it is highly 
likely that willow-associated ECM fungi are crucial for the survival 
of this endangered orchid. The same might apply to two other rare 
orchid species in the Netherlands, Platanthera bifolia and Dactylorhiza 
incarnata ssp. coccinea of which the associated ECM fungi are still 
largely unexplored but which are also mostiy found near Salix 
repens [69]. 

Contribution of Soil Data to Mapping Rare Fungi 

The Netherlands boasts one of the most efficient programs in 
the world for gathering distribution data on macrofungi, where the 
Netherlands Mushroom Mapping Working Group (Werkgroep 
Paddenstoelenkartering Nederland, WPN) has been producing 
distribution maps and fruiting statistics since 1980 with a nation- 
wide network of ca. 600 volunteers that produce appr. 60,000 data 
points yearly [70], [71], [72]. The long-term mapping program of 
the WPN has been crucially important in compiling a near 
complete checklist of macrofungi [73], in highlighting high 
diversity sites [70], and in drawing attention to the serious decline 
of macrofungi in the Netherlands [72]. As a result, an official 
national Red List of fungi was published in 1 996, followed by a 
revision in 2008 [74]. On the other hand, due to inherent 
limitations of the data collecting method, data on the diversity and 
distribution of microfungi and taxa with inconspicuous fruitbodies 
have remained scarce. Building on the methodological advances 
provided by next-generation sequencing techniques, we can still 
greatly improve the exhaustiveness of species inventories by 
providing complementary information, particularly with regard to 
asexual fungi and species with microscopic and/ or inconspicuous 
fruiting bodies, in addition to providing new spatial data points for 
macrofungi already featured in the species lists. 

The massive amount of sequence data we generated provide the 
first kingdom-wide fungal diversity assessment for this coastal 
ecosystem, offering a complementary view on species diversity 
with a particularly large contribution on asexual fungi and species 
with microscopic and/or inconspicuous fruiting bodies that fall 
beyond the scope of the above mapping project. Even for 
macrofungi, despite the obvious spatial limitations of the 
methodology we used, we detected numerous red listed species 
in our samples, often from previously unknown locations. Several 
taxa with scarce (<5) previous sporocarp records, such as Inocybe 
exilis, Pseudobaeospora pyrifera, Russula pascua, Scutellinia vitreola and 
several Tomentella and Sebacina species to name a few, were found in 
multiple locations, indicating that they may be more widespread in 
the Netherlands than previously thought (Figure 4). This may be 
particularly true for Tomentella, a species-rich genus with often 
inconspicuous sporocarps, that is regularly shown to be much 
more diverse and abundant in soil and root samples than in 
sporocarp-based diversity assessments (e.g., [15], [75], [7 6]). We 
also provide valuable data for additional taxa that are considered 
either rare or sporadic, often with "vulnerable" or "endangered" 
red list categories, e.g., Clavicorona taxophila, Clavulinopsis corniculata, 
C. luteoalba, Conocybe subxerophytica, Cortinarius alnetorum, C. casimiri, C. 
cinnamomeus, C. diasemospermus, C. parvannulatus, Entoloma asprellum, E. 
clandestinum, E. longistriatum, E. turbidum, Geastrum striatum, Geoglossum 
fallax, G. umbratile, Geopora arenicola, Hebeloma collariatum, H. leucosarx, 
H. nigellum, Hemimycena tortuosa, Hydnum rufescens, Hygrocybe helobia, H. 
vitellina, Inocybe agardhii, I. decipiens, I. dunensis, I. griseovelata, I.jacobi, I. 



nitidiuscula, I. serotina, I. vulpinella, I. whitei, Lactarius aurantiacus, 
Mycena albidolilacea, Panaeolus papilionaceus, Peziza ampliata, Pleurotus 
eryngii, Russula foetens, Russula puellaris, Squamanita odorata, Tephrocybe 
tylicolor, Thanatephorus cucumeris, Tricholoma lascivum, Tulostoma 
brumale, Xerocomus riperiellus etc. 

Highly similar matching sequences from species not formerly 
reported from the Netherlands represent a special case and 
confirming the presence of some species will require additional 
evidence. Some such examples among ascomycetes are the truffle 
Tuber anniae, the earth-tongue Geoglossum glabrum, the cup fungus 
Geopora cervina, and the false-morel Verpa digital formis. In Basidio- 
mycota, the list of taxa formerly not reported from the Nether- 
lands and found among the highly similar (>98%) matching 
sequences of our OTUs include: Clavaria californica, Cortinarius favrei, 
C. sertipes, Heterobasidion parviporum, Hygrocybe persistens, Inocybe giacomi, 
Mycena simia, Pluteus satur, Russula Integra, R. nana, Sphaerobolus ingoldii, 
and Tubaria minima. Occurrence data regarding species that usually 
occur in arctic-alpine regions, such as Cortinarius favrei and Russula 
nana, are particularly interesting, given the geographic position and 
the lack of mountains in the Netherlands. In the Arctic, these 
species are symbionts of various dwarf willow species (e.g. Salix 
polaris, S. reticulata) [77],[Geml pers. obs.], while in the Nether- 
lands, they appear to grow with Salix repens. Both Cortinarius favrei 
and Russula nana have been found associated with Salix repens in the 
oceanic archipelago north of Scodand [26]. It has to be noted that 
Russula laccata, a species originally described from the Netherlands 
albeit with no current mapping data, is very closely related to R. 
nana and whether or not they are truly distinct species is currendy 
uncertain. Therefore, detailed taxonomic work is needed to test 
whether or not the two taxa are conspecific and whether R. laccata 
refer to the same taxon as our OTU (no. 695), in which case the 
name R. nana will be preserved based on priority. 

Conclusions 

Fungal communities are notoriously difficult to fully character- 
ize for ecological and biodiversity studies and for conservation 
purposes. Even for macrofungi (e.g., such as mushrooms, true and 
false truffles etc.), that have the longest history of diversity studies 
among fungi [1],[70],[78], basic questions about the number of 
species at a given location or differences in species richness among 
vegetation types have generally remained unanswered due to 
taxonomic problems and the scarcity of long-term sporocarp- 
monitoring projects [78]. 

Our work provides a wealth of information on the kingdom- 
wide diversity of fungal communities along the protected coastal 
dune habitats in the Netherlands, extending both the list of species 
likely occurring in the Netherlands as well as providing new data 
points for many rare species. Furthermore, our work highlights 
that fungal communities found in the rhizosphere of S. repens can 
have very different species composition due to edaphic factors and 
that the protection of areas that harbour such varied S. repens beds 
is vital for the conservation of co-habiting fungi, including many 
red-listed species. Future fruitbody collecting projects should 
utilize our data to concentrate sampling and monitoring efforts 
on selected localities in order to confirm species occurrences that 
are currendy only based on soil sequence data. Although more 
detailed, taxon-specific studies will follow, this project already 
provides examples for the potential contribution of high-through- 
put soil sequencing studies to fungal mapping and conservation. 

Data Accessability 

The raw data have been submitted to the European Nucleotide 
Archive (ERP001713). 
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Table SI List of OTUs that showed >97% ITS2 sequence 
similarity to fully identified fungi in coastal Salix repens communities 
used in this study with best identified match (GenBank or UNITE 
accession number, Species Hypothesis, name), similarity (%), 
pairwise alignment length (bp), sequence counts in the sampled 
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